Thermodynamic phase-field model for microstructure with multiple components and 

phases: the possibility of metastable phases 
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A diffuse- interface model for microstructure with an arbitrary number of components and phases 
was developed from basic thermodynamic and kinetic principles and formalized within a variational 
framework. The model includes a composition gradient energy to capture solute trapping, and is 
therefore suited for studying phenomena where the width of the interface plays an important role. 
Derivation of the inhomogeneous free energy functional from a Taylor expansion of homogeneous free 
energy reveals how the interfacial properties of each component and phase may be specified under 
a mass constraint. A diffusion potential for components was defined away from the dilute solution 
limit, and a multi-obstacle barrier function was used to constrain phase fractions. The model was 
used to simulate solidification via nucleation, premelting at phase boundaries and triple junctions, 
the intrinsic instability of small particles, and solutal melting resulting from differing diffusivities 
in solid and liquid. The shape of metastable free energy surfaces is found to play an important role 
in microstructure evolution and may explain why some systems premelt at phase boundaries and 
phase triple junctions while others do not. 
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I. INTRODUCTION 

Developing an understanding of microstructure forma- 
tion in multiphase, multicomponent systems is a chal- 
lenge for industrial development of advanced alloys, yet 
interesting from a philosophical perspective due to the 
formation of complex patterns for which no theory exists 
[1]. Compared to two-phase binary systems, multiphase 
and multicomponent systems have additional degrees of 
freedom that introduce inherent complexity. Multiphase 
systems have the ability to form phase triple junctions, 
transient phases [2 , and metastable phases at grain 
boundaries or triple junctions (i.e. interfacial premelt- 
ing) |3114j. 




FIG. 1: A hypothetical ternary phase diagram. A and 

B denote the bulk composition of two phases in 

equilibrium. Diffuse interfaces are curves that connect 

A and B. 

Additionally in multicomponent systems, adsorption 
of components to interfaces or triple junctions is possi- 
ble. Consider figure [l] which illustrates a hypothetical 
ternary phase diagram. A and B denote the bulk compo- 



sition of two phases in equilibrium, and a diffuse interface 
between A and B corresponds to a curve connecting the 
points. At equilibrium there can be only one curve. In 
binary systems, the interfacial profile is constrained to lie 
on the dotted line that represents a linear combination 
of the components. In a ternary system however, the ad- 
ditional component permits complex pathways that de- 
pend on the free energies of the phases, the presence 
of metastable phases, and the energy-minimizing path 
through composition space. 

In this work we develop a multicomponent, multiphase 
model that treats the diffuse interface in a thermody- 
namically consistent way, allowing us to investigate pre- 
melting and the effects of metastable phases on inter- 
facial composition in multiphase systems. The model 
includes a (Vc)^ in the free energy functional as a nat- 
ural way to model the correct amount of solute trap- 
ping [5 . We present a careful derivation of the model in 
order to cast the multiphase, multicomponent problem 
within a thermodynamic framework, and derive compo- 
nent diffusion equations that obey the Gibbs-Duhem and 
Nernst-Einstein relations. Subtle differences are clarified 
between the chemical potential and the diffusional poten- 
tial. These differences become important in multicompo- 
nent systems. Specifically, the driving force for diffusion 
is sometimes defined as |^ [6,-^ , but this definition is in- 
correct for a multicomponent system that obeys a mole 
fraction constraint. 

Phase-field has emerged as an important method for 
modeling microstructure evolution because of its abil- 
ity to simulate complex geometries while incorporating 
thermodynamic and kinetic data. A phase-field model 
assumes that interfaces in microstructure are diffuse at 
the nanoscale and can be represented by one or more 
smoothly varying order parameters, eliminating the need 
to explicitly track boundaries. Nonlinear diffusion and 
curvature-driven physics are incorporated, and creation. 



destruction, and merging of interfaces are handled implic- 
itly. A substantial amount of literature has been written 
on phase- field models and is summarized in recent review 
papers [TQUT4]. 

Phase-field models may be classified into two categories 
based on their philosophical treatment of a diffuse inter- 
face. In the approach pioneered by Cahn and Hilliard 
[15j , the interface is a coarse-graining of the underlying 
atomistic representation. The width of the interface in 
the model is identical to its physical width, which may 
be as small as a few nanometers. This approach pro- 
duces a thermodynamically consistent description of the 
interface, but makes simulation of realistically sized mi- 
crostructures problematic. Microstructural features are 
often on the order of micrometers or larger, several or- 
ders of magnitude larger than the interfacial width. This 
presents a computational challenge for simulating large 
microstructure. In principle, the problem could be ad- 
dressed with faster computers and improved numerical 
algorithms. 

The second approach uses phase-field as a modeling 
tool for solving the underlying free boundary problem 
without explicitly tracking boundaries. It is numeri- 
cally advantageous to allow the computational interfacial 
width W to exceed the physical width, but doing so in- 
troduces error that scales with W [16]. The "thin inter- 
face" approach [17, 18 was an important development 
in this regard. With the appropriate choice of model 
parameters, thin interface models converge to the Gibbs- 
Thomson relation in the limit where the interfacial width 
is much smaller than a typical pattern size of the system. 
As a result, convergence is on the order of W^. Ex- 
cessive solute trapping occurs when the numerical width 
becomes large, but has been remedied with anti-trapping 
currents [19 . Notably, this approach has produced sim- 
ulations of dendrites that are quantitatively comparable 
to experiment. 

Alloy phase- field models have been developed following 
both philosophies and will be briefly reviewed. Wheeler, 
Boettinger, and McFadden (WBM) [H [20l ET treated 
the interface in a thermodynamic way but were limited 
to binary systems with two free energy curves due to 
fundamental model difficulties. Steinbach et al. ^22i i23| 
prompted development of a series of models for multi- 
component and multiphase systems that have produced 
quantitative simulations on experimental length scales 
[24] . However these multiphase models are not appropri- 
ate for studying phenomena where the interfacial width 
plays an important role, such as solute trapping, inter- 
face premelting, nucleation, or the appearance of tran- 
sient phases. 



A. The Wheeler-Boettinger-McFadden model 

The Wheeler-Boettinger-McFadden (WBM) model [5, 
[201 121] begins with the Cahn-Hilliard free energy func- 
tional for a binary system [TB^ (see section II A[) and in- 



troduces a non-conserved order parameter (f) to indicate 
which regions of the system are solid (^ = 1) and which 
are liquid {(j) = 0). At an interface between liquid and 
solid, both (j) and c vary smoothly from one phase to the 
other. The free energy functional for the system is: 



F[c,4>] 



fo{4>, c, T) + -ec(Vc)2 + -e^{V4>f dV 

(1) 
The homogeneous free energy density fo{(p^c^T) pro- 
motes phase separation in the absence of interfacial en- 
ergies, and Cc and e^ are the composition and phase gra- 
dient energies, respectively. Phase and composition gra- 
dients overlap at equilibrium to form an interface, and 
the gradient squared terms smooth the interface and in- 
troduce interfacial energy. The (Vc)^ term was omitted 
from the original model for computational convenience 
[20 , but was later found to be necessary for modeling 
solute drag during rapid solidification [5 . 

The WBM approach models a diffuse interface as an 
interpolation between phases where the composition of 
phases at the interface are equal. An interpolating func- 
tion p{(j)) is used to connect the homogeneous free energy 
densities of the phases: 

fo{4>, c, T) = p{<P)fHg{c, T) + (1 - p{<P))fsoi{c, T) (2) 

Interpolation between two free energy curves is illus- 
trated in figure |2b[ p{(j)) has a minimum at </> = and 
(j) = 1 and provides a barrier for transition from one phase 
to the other. It lacks a physical basis and is generally cho- 
sen for numerical convenience. The WBM model requires 
that p(0) approach = and = 1 with zero slope in 
order to prevent the appearance of negative phase frac- 
tions. 

Because there is no natural extension of p{(j)) to han- 
dle an arbitrary number of free energy curves, there have 
been few attempts to develop a multiphase model follow- 
ing the WBM approach. Folch and Plapp [16 derived a 
thin-interface model that included an interpolating func- 
tion for three curves, and Nestler et al. [25 developed 
a WBM-like nonisothermal multiphase, multicomponent 
model. However, neither work included a composition 
gradient energy or applied the correct thermodynamic 
constraints to the component diffusion equations. These 
issues are addressed in sections [TTl and ITITl 



B. The Access multiphase model 

Steinbach and co-workers developed the so-called Ac- 
cess multiphase model, the first phase-field model capa- 
ble of simulating the interaction of an arbitrary number 
of phases [22, 23 . The original model did not include 
solute diffusion and considered pairwise interactions be- 
tween phases using double well interpolation functions 
and Allen-Cahn dynamics. Modeling the dynamics of a 
multiphase system as the sum of pair-wise interactions 
was problematic, violating Young's Law at phase triple 



junctions, and was fixed with the introduction of inter- 
face fields [23] , 

Tiaden et al. [26] added single-component solute dif- 
fusion to the Access multiphase model. The interfa- 
cial region was modeled as a blend of phases, each with 
a phase fraction (po^ and unique composition c^, con- 
strained so that the concentration of the system was 
c(x, t) = ^^ (j^aCa- The diffusing species was partitioned 
amongst the different phases, and Fickian diffusion equa- 
tions were solved in each phase. The diffusion equations 
were coupled to phase evolution equations driven by a 
difference in free energy between phases which was de- 
termined from a local linearization of the phase diagram. 

The dilute solution limitation of the Tiaden model was 
removed in an extension by Kim et al. [27 for single- 
component diffusion with the use of an interpolating 
function to link the free energy curves. Kim also intro- 
duced a more sophisticated condition of equal chemical 
potential to determine how to distribute solute amongst 
the phases at a diffuse interface. 

Grafe et al. ^28 developed the first multicomponent 
extension of the Tiaden model. The driving force for dif- 
fusion was Vcf , the concentration gradient of component 
i in phase a, which is a dilute solution approximation. 
Solute distribution was calculated with partition coeffi- 
cients from Thermo- Gale. 

Eiken et al. [29] developed a multicomponent exten- 
sion of the Tiaden model which removed the dilute solu- 
tion limitation and allowed for easier inclusion of thermo- 
dynamic data. Jl^ = ^^ was chosen as the driving force 
for diffusion, although Jlf is the slope of the free energy 
curves and not the chemical potential /if = |^ (see 
section III A). A very computationally expensive quasi- 



equilibrium calculation was necessary at each timestep 
to relieve the dilute solution approximation. 

A two-phase multicomponent model with an antitrap- 
ping current was presented by Kim [30 , but to date an 
antitrapping current has not been included in a model 
with both an arbitrary number of phases and compo- 
nents. 




(a) Molar free energy for two phases, a and /3. 




(b) An interpolating function is used to smoothly 
connect the molar free energy curves. 

FIG. 2: The WBM model assumes that the energy of 

interfacial compositions is a weighting of the dashed 

regions of the free energy curves, while Access models 

assume that the energies of interfacial compositions lie 

on the common tangent line. 



C. Graphical interpretation 

The fundamental difference between the WBM and Ac- 
cess approaches is illustrated in figure |2a| where free en- 
ergy curves and the common tangent construction for 
two phases a and /3 are drawn. A diffuse interface must 
include compositions between the equilibrium concentra- 
tions Cgg and cfg, but the energy of these intermediate 
compositions is somewhat ambiguous. 

The WBM model assumes that each phase at an inter- 
face has the same composition. At equilibrium, the free 
energy of these interfacial points is then a weighted av- 
erage of the dashed portions of the free energy curves in 



figure |2a| The dotted line in figure [2a| indicates a poten- 
tial path when a barrier in (j) is added, and corresponds 



to the dotted path lying on the free energy surface in 
figure 2b When the system is not at equilibrium, in- 



terfacial compositions may lie anywhere on the surface 
of figure |2b[ A/ denotes energy at an interfacial point 
relative to a composite blend of a and /3 at equilibrium. 
The gray shaded region is A/ integrated across an inter- 
face, and represents the interfacial energy contribution 
from including intermediate compositions at the inter- 
face. The contribution of the shaded area increases for 
wider interfaces because more material with energy above 
the common tangent construction must be introduced. 

The Access approach for modeling interfaces is to as- 
sume each phase has its own unique composition that 
cannot be measured experimentally but which evolves 
toward its equilibrium concentration. Interpolation be- 
tween phases at their equilibrium concentration produces 



intermediate compositions with energies that he on the 
common tangent hne with a barrier only in ^, prohibit- 
ing the appearance of metastable phases which he above 
the common tangent. 



The gray hne in figure 2b illustrates an Access inter- 
face that connects the common tangent points of the free 
energy curves. Because the barrier in (j) is the only con- 
tributor to A/, widening the interface for computational 
convenience does not introduce more material with en- 
ergy above the common tangent. 



II. THE MULTIPHASE FREE ENERGY 
FUNCTIONAL 

A. Free energy of a binary system 

In an infiuential paper that laid the foundation for 
phase-field modeling, Cahn and Hilliard derived an ex- 
pression for the free energy of an inhomogeneous binary 
system [15 . Their approach was to assume that the free 
energy of an infinitesimal volume in a nonuniform system 
depends both on its composition and the composition of 
its nearby environment. Total free energy cannot depend 
solely on local composition because different spatial con- 
figurations with the same volume fraction are not en- 
ergetically equivalent; a heterogeneous system has more 
interfacial area and will have a higher energy. 

Starting with the homogeneous free energy density for 
a binary system /o(c), they performed a Taylor expansion 
in terms of the derivatives of composition to approximate 
/(c, Vc, V^c, . . .). Morris [31[ provides justification for 
excluding terms linear in |Vc|. For isotropic or cubic 
symmetry {df/d\/c)o = 0, and the free energy simplifies 
to an equation with constant coefficients and even powers 
of Vc: 

/ = /0(c)+/.iV'c+i/.2(Vc)2 + iA.3(V'c)2+/.4V^C+... 

(3) 
Cahn and Hilliard then argued that the derivative terms 
with even powers V^c, V^c, V^c, etc. should vanish. Be- 
cause of the assumption that the free energy density is 
inffuenced only by concentration within a small neighbor- 
hood, it is reasonable to truncate the expansion. Keep- 
ing only terms up to second-order produces the Cahn- 
Hilliard free energy functional: 



F[c] 



Ik 
Jv 



(c) + k(Vc)2 dV 



(4) 



where z^: is a gradient energy coefficient that penalizes the 
formation of sharp interfaces. 



B. Free energy of a multicomponent system 

The approach of Cahn and Hilliard is now applied to 
a system with an arbitrary number of components. A 



system with M components has M — 1 independent mole 
fractions that obey the following constraint: 



M 



(5) 



The inhomogeneous free energy becomes a function of 
each independent component as well as their derivatives: 

/(ci,C2,...,Vci,Vc2,...,V2ci,V2c2,...). The Taylor 
expansion [32 of the multicomponent / about a homo- 
geneous point /o = /(ci, C2, . . . , 0, 0, . . . , 0, 0, . . .) is: 

/(Ci, C2, . . . , Vci, VC2, . . . , V^Ci, V^C2, . . .) 
= /o(ci,C2,...) 



(a^)/'- + (3^)/'= 



0>V2 



1 



,2 



V^ci 



df 



av2 



V^ca 



C2/0 






(VC2)2 + ... 



(6) 



For simplicity only terms for two components up to 
second-order have been written out because higher order 
terms will be excluded. Once again only isotropic and 
cubic symmetry of / is considered, allowing the tensors 
to be replaced with constants, and terms in Vc, as well 
as even derivatives, are excluded. The Taylor expansion 
of f in terms of the M — 1 independent components and 
their derivatives is: 



M-l 



M-1 



/ = /o(ci,C2,...)+ ^ ^l^i{'^Ci)^^Yl X] ^ijVCi-Vc^- 

(7) 
This equation, which has been previously reported in lit- 
erature [6^, T", '33^, may be simplified by combining the 
summation terms: 



M-l 



/ = /o(ci, C2, . . .) + Yl o'^^i^^^ • ^^. 



(8) 



*,J=1 



where K,ij is a symmetric matrix of gradient energy coeffi- 



cients discussed in section II F The free energy functional 
for a multicomponent system is: 



r M-l 

F[{c}]= /o({c})+^-/..,Vq-Vc,) 



dV (9) 



where {c} denotes a set of M — 1 independent mole frac- 
tions. 



C. Definition of a phase 

A phase is a region of a microstructure with homoge- 
neous properties that is physicahy distinct from other re- 
gions of the system, excluding geometric transformations 
that map one region onto another. Phases in microstruc- 
ture commonly differ in composition and/or crystal struc- 
ture, although many other physical differences are pos- 
sible. The volume fraction of phases in equilibrium is 
predicted from thermodynamics, but phase itself is not a 
thermodynamic state variable; phase is a labeling device 
that identifies a unique thermodynamic state function. 

Each phase a is assigned a phase fraction (po^ that varies 
between and 1. (j)^ = ^ designates areas where no a- 
phase is present, and (j)a = I corresponds to single-phase 
regions of a. For a system with N phases, the phase 
fractions obey the following constraint: 



N 



Y.K = l 



(10) 



a=l 



Microstructure (excluding grain boundaries, defects, 
etc.) is composed of single phase regions separated by 
interfaces, and only at interfaces are more than one (j) 
non-zero. The interface between two phases is assumed 
to consist of a thin layer across which the physical prop- 
erties vary continuously from one phase to the other, 
and a diffuse interface at equilibrium represents a bal- 
ance between free energy curves, composition gradients, 
and phase gradients. 

Because the thermodynamic potential of a multiphase 
system is equal to the summation of potentials over all 
phases, a linear weighting of the free energy densities by 
phase fractions is used to represent the homogeneous free 
energy of a multiphase system: 

N 
/0({C}, {(/.I, 02, . . . cPn}) = Y, ^cUic}) (11) 



This form reduces to /a({c}) when only the a-phase is 
present, yet can be constructed for an arbitrary number 
of phases. 
Eq. 



11 does not provide an energy barrier for diffu- 



sionless phase transformations, and without modification 
does not correctly describe phase transitions of pure com- 
ponents. A simple barrier between a and (3 of the form 
^a/3^a^/3 IS Suggested, although the multi-obstacle bar- 
rier introduced in section |IV A| permits any function to 
be used as a barrier. The homogeneous free energy be- 
comes: 

N 
fo{{c},W},T) = Y,^afa{{c},T) + J2Wap^e.'Pp 

(12) 
Wai3 > captures the mean field interaction between 
phases and is analogous to a positive enthalpy of mixing 
for phases. 



D. Free energy of a multiphase, multicomponent 
system 

For an A^-phase, M-component system, 

/({c}.{^}.{Vc},{V^},...) is a function of M - 1 
independent mole fractions and A' — 1 independent 
phase fractions. Because / can describe nonequilibrium 
systems where metastable phases are present, it is not 
necessary that the Gibbs phase rule be obeyed. 

Once again isotropic and cubic symmetry of the free 
energy is considered, terms in Vc and even derivatives of 
c are excluded from the Taylor expansion, and only terms 
up to second-order are kept. The full expansion about 
the homogeneous free energy /o({c}, {(/>}, {0}, {0}, . . .) is 
not algebraically difficult but has many terms and is not 
explicitly written out here. It is analogous to Eq. [9] but 
with two additional sets of terms. One set couples phase 
gradients with gradient energy coefficients X^fd- The sec- 
ond set couples composition gradients and phase gradi- 
ents: 



dy 



dVcidV4)a 



VCi ■ V<^a 



(13) 



The coefficients of these terms form an M x A^ matrix 
(_ia, and introduce an additional energy penalty for over- 
lapping phase and concentration gradients. The total 
gradient energy contribution for a multicomponent, mul- 
tiphase system may be written in compact form as: 



[Vc V^] 






a/3 



Vc 

V<^ 



(14) 



For simplicity we assume that ^ia = in this work. 
The multiphase, multicomponent free energy functional 
then becomes: 



r N-1 

-^^ L a, (3=1 



M-1 



(15) 



E ^/^ijVQ-V^ 



*,j=i 



dV 



The free energy curves are the driving force for phase sep- 
aration, and the gradient energy coefficients Xa(3 and hZij 
penalize gradients that develop, creating a surface energy 
at phase boundaries. f<iij penalizes phases for differing in 
composition, and Xa(3 introduces additional energy not 
captured by the composition gradients at phase bound- 
aries. This energy derives from some physical difference 
between the phases other than composition. Eq. [T5J 
which is the central equation of focus in this work, is a 
first order approximation of the free energy of a system 
with an inhomogeneous distribution of phases and com- 
ponents. It reduces to the Cahn-Hilliard equation for a 
two-component system. 



E. Surface energy 

The surface energy of a diffuse interface is the excess 
grand canonical potential, the difference between the free 
energy functional and the minimized free energy the sys- 
tem would have if the properties of the phases were con- 
tinuous: 



(7 = minF[{c}, {(/)}] 






Mi^i 



(16) 



The first term is the minimum of F found by application 
of the Euler-Lagrange equation, and ^^ /i^Q is the ho- 
mogeneous free energy /if is the chemical potential of 
component i at equilibrium and is found by computing 
the tangent plane to the free energy surfaces. 

Surface energy in this model has two contributions. 
One contribution comes from phase and composition gra- 
dients which are present at the interface, and the other 
results from composition deviating from its equilibrium 
value at the interface as illustrated in figure l2al 



F. Interpretation of the gradient energy matrices 



Although Eq. 15 is a function only of the independent 
gradients, it is necessary to specify the properties of the 
dependent component and phase as wellQ For an TV- 
phase, M-component system, the phase gradient energy 
matrix A has A^ — 1 rows and columns and the component 
gradient energy matrix K: has M — 1 rows and columns. 
The gradient energy coefficients coupling the implicitly 
defined N^^ phase (and M^^ component) are not explic- 
itly defined in A and />:, but are instead distributed across 
all of the coefficients. A and Hi are dense versions of larger 
matrices, A and K, that have a direct physical interpre- 
tation. The complete coupling of all gradients can be 
written in matrix form: 



[V(/)i V(/)2 ••• V(t)N] 



All Ai2 ••• AiatI rV(^i 

A21 A22 • • • A2Ar V^2 



Aati AAr2 ••• AatatJ [vc^at 

(17) 
The coefficients of A specify an energy penalty for ev- 
ery possible pair of overlapping gradients. An analo- 
gous M X M matrix K contains composition gradient 
energy coefficients Kij that penalize overlapping compo- 
sition gradients. 

If the phase conservation constraint V^at = — (V(/)i + 
V(/)2 + . . . + \/(j)N-i) obtained from Eq. 10 is substituted 



into Eq. 17 and the matrix multiplication is performed, 
an expression representing the gradient energy in terms 
of the A" — 1 phase gradients is obtained. The coefficients 
in this expression are related to the Xa(3 that form the 
matrix AJ^ Because of the dependence of the A^^^ phase 
on all other phases, elimination of the A'^^ row and col- 
umn of A distributes the gradient energy coefficients for 
the A^^^ phase across all coefficients of A. Thus A will 
generally be a fully dense matrix. 

The physical basis for A and K requires that they be 
symmetric positive definite matrices. A and K must be 
positive-definite because if they had negative eigenvalues, 
there would be a coupling of gradients (in the direction 
of the corresponding eigenvector) for which an increas- 
ingly sharp interface lowers the free energy of the sys- 
tem, producing a physically impossible negative surface 
energy and rendering the evolution equations unstable. 
The simplification to reach Eq. [8] and Eq. 15 also reveals 
that A and k are symmetric. 



G. Gradient energy coefficient selection 

The free energy functions and gradient energy matrices 
are coupled by Eq. [15] in a way that makes fitting gradi- 
ent energies to experimental systems potentially cumber- 
some. At equilibrium, the surface energy of an interface is 
fixed once the free energy densities and interfacial widths 
are specified. Thus in principle the gradient energy co- 
efficients could be obtained by measuring the width of 
both the composition and phase variations at all possible 
equilibrium interfaces. For a system with A" phases and 
M components, there are potentially (^) = ^ ( A'^ — A^) 
unique phase interface widths, and (^) = |(M^ — M) 
unique compositional interface widths. The number of 
unique widths correspond exactly to the number of up- 
per diagonal coefficients in k. and A. 

However, the shape of free energy functions may pre- 
clude the formation of many possible interfaces, making 
it impossible to determine gradient energy coefficients 
from equilibrium observations. In this case some gradi- 
ent energies takes on a non-equilibrium role, and it may 
be possible to fit the coefficients to equilibrium interface 
widths by assuming that some of the gradient energies 
are zero. In the general case that all gradient energies in 
are non-zero, determination of the coefficients is 



Eq. 14 



nontrivial. The number of unique coefficients is signifi- 
cantly larger than the number of equilibrium observables. 
A series of ab initio calculations would be necessary to 
determine the coefficients. For each coefficient, the in- 
crease in energy when a homogeneous system is forced to 
incorporate a gradient must be calculated. 



^ We found that ignoring the dependent phase and component 
produced an unexpected asymmetry in interfacial compositions 
that was visible in composition maps like those in figure [6c] A 
correct treatment of the gradient energy matrices removed the 
asymmetry. 



^ The diagonal terms Xaa are the coefficients of the squared terms, 
and the off diagonal terms X^p are equal to the coefficients of 



the cross terms multiplied by ^. 



III. COMPONENT EVOLUTION 

Component evolution equations are derived here for a 
non-ideal ternary system. Extension to a different num- 
ber of components follows the same approach but is alge- 
braically tedious. Parts of this derivation are drawn from 
work by Nauman and Balsara [34] and Nauman and He 



The thermodynamic condition defining equilibrium in 
phase-separating systems is the elimination of all chem- 
ical potential gradients. Fickian diffusion with Vc as 
the driving force applies only to the special case of an 
ideal system where there is no enthalpy of mixing. Sys- 
tems which undergo phase separation exhibit "uphill dif- 
fusion" , and Fickian diffusion does not hold. 

To derive component evolution equations for a system 
characterized by a free energy functional, it is necessary 
to begin with the generalized form of Pick's first law: 



Ji = -MiVfii 



(18) 



Ji is the flux of component z, Mi is its mobility, and (li is 
its diffusion potential. In principle Ji might also depend 
on the diffusion potential gradients of other components 
besides i, but this is not considered here. Mi is related 
to the diffusivity Di by the Nernst-Einstein relation: 



M, 



DiC 



RT 



(19) 



If Di depends weakly on composition, q will be the lead- 
ing term in the mobility expression. It is important that 
mobility depend on q for conserved quantities. If it did 
not, it would be possible to have a flux of a component 
without any of that component being present initially. 



the standard definition of chemical potential: 

'9(minF)' 



Mi 



dG 

m 



T,P,Nj^, 



dNi 



(21) 



T,V,Nj^, 



The inhomogeneous chemical potential is defined away 
from equilibrium and approaches the classical chemical 
potential as equilibrium is approached. At equilibrium 
F is minimized, fii is no longer a function of position, 
and fli = jii. 

The fundamental relation for the ternary free energy 
functional F at constant temperature and pressure can 
be written as: 



F 



F 

N 



-PV + ftiCi + /i2C2 + ftsCs 



(22) 



where F is a molar quantity, V is molar volume, and 
N = Ni -\- N2 -\- Ns is the total number of moles in the 
system. It can be shown by standard thermodynamic 
arguments that the fli obey a generalized Gibbs-Duhem 
relation at constant temperature: 



y^Cidfii = VdP 



(23) 



Application of the mole fraction constraint (Eq. p| to 
Eq. [22] to eliminate C3 reveals that the variational deriva- 
tives of F with respect to q are related to differences in 
inhomogeneous chemical potentials: 






Ml -Ms 



(24a) 



T,y,c2 



A. Generalized diffusion potential 

The free energy functional F is a non-equilibrium gen- 
eralization of Helmholtz free energy that includes contri- 
butions from concentration and phase gradients. At equi- 
librium the functional is equal to the equilibrium free en- 
ergy. To describe kinetic evolution in a non-equilibrium 
system, it is necessary to define a potential that ap- 
proaches the chemical potential at equilibrium. Since F 
is a functional, the functional derivative defines an inho- 
mogeneous (or variational) chemical potential field that 
becomes uniform at equilibrium: 



Mi 



5F 



SNi{x) 



(20) 



T,V,N,^i 



Ni (x) is the number of moles of component z as a function 
of position|j Hat notation indicates that the inhomoge- 
neous chemical potential fli is a different quantity from 



^ Throughout the rest of this paper, variational derivatives will be 
written with the assumption that the function in the denomina- 
tor depends on x. 



6F\ 



M2 -M3 



(24b) 



T,V,ci 



These quantities are diffusion potentials, and may be 
interpreted as the energy change upon adding a small 
amount of q while simultaneously removing a small 
amount of C3. Thus the equilibrium condition of con- 
stant inhomogeneous chemical potential is equivalent to 
constant diffusion potential for a system with a mass con- 
straint. 



B. Evolution equations 

The derivation of evolution equations begins with the 
observation that when individual chemical potentials are 
defined, their gradients are related by the Gibbs-Duhem 
equation. If local thermodynamic equilibrium is as- 
sumed, the Gibbs-Duhem equation relates gradients in 
chemical potential V/ii instead of changes in chemical 
potential dfli. Local equilibrium implies that global in- 
tensive parameters vary so slowly that small neighbor- 
hoods around a point can be considered at equilibrium. 
Furthermore for solids and liquids, VdP is generally very 



small and can be neglected in Eq. 23 for simplicity. The 
Gibbs-Duhem relation for an inhomogeneous ternary sys- 
tem then becomes: 



Ci V/ii + C2V/i2 + CsVjls 







(25) 



The mole fraction constraint (Eq. Isl) is used to eliminate 
C3, and the equation is rearranged to put V/ii on the left 
hand side: 

V/ii = (V/ii - V/is) - ci(V/ii - V/ia) - C2(V/i2 - V/ia) 

= (1 - ci)V(/ii - As) - C2V(/i2 - As) 

(26) 



The variational derivatives (Eq. 24) can now be substi- 
tuted in place of the chemical potential differences: 

SF 5F 

V/ii = (1 - ci)V- C2V—- (27a) 

OCi 0C2 

A similar procedure is used to find V (12'- 

6F 5F 

V/i2 = (1 - C2)V-^ - ciV-^ (27b) 

OC2 OCi 



The dynamics of component diffusion is governed by a 
mass conservation law: 



dcj 

dt 



-V-Ji 



(28) 



Substitution of Eq. Tsj T9J and 27 produces component 
diffusion equations for a ternary system: 



dt 



Die ( SF 5F 



(29a) 



do 



.v.r^fn 



SF SF\\ 

The variational derivative ^ is found by applying the 
Euler-Lagrange equation to the free energy functional 
(Eq. [15]): 



SF SF ^^ dfa Y^ ^2 fnn\ 



IV. PHASE EVOLUTION 



where Ta is a kinetic coefficient associated with how 
quickly the a-phase can transform to another phase at 



constant composition. 



SF 



is found by applying the 



Euler-Lagrange equation to the multiphase multicompo- 
nent free energy functional (Eq. IlS]): 



SF dfo 



N-1 
/3=1 



(32) 



with /o defined in Eq. [12] The implicitly defined phase 
fraction (J)n is a function of the other phase fractions 
such that -^^ = — 1. Thus the driving force for phase 
separation becomes fa — In^ where /n is the free energy 
density of the implicitly defined N^^ phase. 



A. A barrier function for phase fractions 

The definition of the phase fraction as a positive quan- 
tity less than or equal to 1 imposes a constraint on the 
phase evolution equations which was not included in their 
derivation. In fact, negative phase fractions would be en- 
ergetically favorable if they had physical meaningj^ Con- 
sider a single component system with a high energy phase 
a and a low energy phase /3. Converting a to P decreases 
free energy by //3 — fa- If negative phase fractions are 
not prohibited, there is an arbitrage where simultane- 
ously producing more (3 and negative a lowers the free 
energy without violating the phase fraction constraint. 
The problem is that the global energy minima are un- 
bounded in (/). Phase-field models typically address this 
issue by constructing /(c, </>) so that it has minima at 
= and 0=1 and penalizes < and > 1. How- 
ever, constructing such an interpolation for an arbitrary 
number of phases is problematic. 

Barrier methods are often applied to minimization 
problems subject to inequality constraints. Constrained 
optimization consists of minimizing the original poten- 
tial plus the barrier functions representing the inequality 
constraints. Logarithms^ and ^ functions are commonly 
used barriers, but are not ideal candidates for phase frac- 
tions which spend a lot of time in the vicinity of = 
where the barriers are undefined. Single phase-regions in 
a multiphase system would be unstable for instance, as 
would any evolution directed along the boundary of the 
feasible region, corresponding to a phase transition. 

A multi-obstacle barrier [25l |38l [39] is used here to 
constrain phase fractions. The barrier is zero for permis- 



Phase fractions are not coupled by thermodynamic re- 
lationships and are not conserved quantities since phases 
are created and destroyed during phase transitions. Thus 
phase evolution follows Allen-Cahn dynamics |37|: 



SF 



dt 



(31) 



^ In systems where borrowing is allowed, negative percentages have 
meaning. Consider financial leveraging - taking out a loan to 
make an investment. It could be profitable to say, invest 150% 
of your income by taking out a -50% loan, if you expect the 
return on the investment to be higher than the interest due on 
the loan. 

^ The cln(c) terms in the ideal entropy of mixing are a barrier 
function for components that has a thermodynamic justification. 



(^2=1 








FIG. 3: The multi-obstacle projection for a three-phase 

system. Black points indicate initial locations in 

phase-space, and gray arrows a possible trajectory in 

the absence of constraints. Constrained evolution 

proceeds along the bent black arrows. 



Generalization of the projection procedure for an N- 
phase system involves fixing violations and then recur- 
sively projecting the system to lower dimensions to fix 
additional violations. The implementation of this recur- 
sive procedure is presented in algorithm [l] 

Algorithm 1 multiObstacle({^i, . . . ^at-i}) 

for 0^ = 01 . . .0Ar_i do 
if (/)i < then 

0^ = O 

end if 
end for 

(/)Ar ^ 1 - YliJl ^i 

if (j)N < then 

for (l)i = (1)1 ... (l)N-2 do 



end for 

multiObstacle({(/)i, . 

end if 



• 0Ar-2}) 



sible phase fractions, and infinite otherwise. The multi- 
obstacle barrier is a generalization of the double obstacle 
barrier used in Access models. The double-obstacle po- 
tential was studied by Blowey and Elliott [38 and found 
to be consistent with curvature dependent phase bound- 
ary motion in two-phase systems. An algorithmic imple- 
mentation of the barrier is presented here for a system of 



A^ order parameters that obey a constraint like Eq. 10 



The phases in an A/'-phase system form the vertices of 
an A^-simplex, and the feasible set of phase fractions lie 
on or within this simplex. Enforcing that all N phase 
fractions remain positive is enough to insure that all N 
phase fractions will also be less than 1 because of the 
phase fraction constraint (Eq. 10). The multi-obstacle 



barrier is implemented by projecting a vector of phase 
fractions back onto the surface of the simplex when one 
or more phase fractions become negative as a result of 
advancing the evolution equations. For a 2-phase system 
there is only one independent phase fraction, and the 
projection is trivial. If ^ < set (/) = 0, but if (/) > 1 set 

4>=i. 

Figure [3] offers a geometric description of the projection 
for a three-phase system. Orthogonal axes are drawn to 
represent the two independent phase fractions ^i and 
02, and each coordinate in the graph corresponds to a 
unique point in phase space. ^3 = 1 at the origin, and 
03 = corresponds to the dashed line connecting 0i = 1 
and 02 = 1- The constraint that all three be positive 
restricts the feasible region to the triangle with vertices at 
the origin, 0i = 1, and 02 = 1. Any point outside of this 
triangle is non-physical and is given an infinite energy 
penalty by setting the offending phase fractions to zero. 
In the case that 03 becomes negative, the projection is 
accomplished by moving in the (—1,-1) direction until 
03 becomes zero. 



V. RESULTS 

Experiments have repeatedly shown that liquids can 
often be supercooled before they solidify, but solids can 
almost never be superheated [^ |40J |41J • Solids often be- 
gin to melt below the bulk melting temperature, with 
liquid appearing first at triple junctions and then at 
grain boundaries [3 . Explanations have included the ob- 
servation that grain boundaries and triple junctions are 
high energy sites that are less thermally stable than the 
bulk [42 , that free surfaces may premelt due to atomic 
thermal vibrations [40], and that premelting my result 
from a structural transition [43] . The results in this sec- 
tion demonstrate premelting in nanostructures due to the 
shape and position of met ast able free energy surfaces, 
possibly explaining why some experimental systems form 
stable liquid films at phase boundaries while others do 
not. 

A four-phase ternary eutectic free energy landscape 
was developed from a ternary regular solution model of 
the form: 



/(Cl, C2) =^12ClC2 + ^13ClC3 + ^23C2C3 

^RT (ci In(ci) + C2 ln(c2) + C3 ln(c3)) 



(33) 



where C3 = 1 — ci — C2, and Vt determines the enthalpic 
contribution to free energy. The function has a minimum 
at (ci,C2,C3) = (^, ^, |), and phases with different equi- 
librium compositions were created by translating Eq. 33 



/i(ci,C2) = /(ci + 1/3 - .9,C2 + 1/3 - .05) 
/2(ci,C2) = /(ci + l/3-.05,C2 + l/3-.9) 
/3(ci, C2) = /(ci + 1/3 - .05, C2 + 1/3 - .05) 
/4(ci,C2,T) = /(ci,C2) + .5266 + A/™ 



(34) 
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These four surfaces are plotted in figure |4j The system is 
a ternary eutectic in the sense that a silver liquid phase 
appears in the center of the phase diagram above the 
melting point, and upon cooling, separates into three 
solid phases, each with a limited amount of solubility. 
The energy of the liquid surface minimum relative to 
the other surfaces is specified by A/^, the free energy 
change upon melting of the solid|j The liquid surface 
is calibrated so that the minima of all four surfaces lie 
on a common tangent when Afm = 0, corresponding to 
T — T 



Although the liquid surface in [4a| lies above the con- 
vex hull of the solid surfaces, there is a small region of 
composition space at the center of the ternary triangle 
where the liquid surface is lower in energy than any of 
the solid surfaces at the same composition. Liquid in 
this region of composition is metastable with respect to 
phase-separation into the three solid phases. The simula- 
tions that follow examine the effect of such a metastable 
region on microstructure evolution. 



The evolution equations (Eq. 29 and 31) were nondi- 
mensionalized as follows: 



Di = D 



T 

tRT 



V 



= n. 



A 



■a^ 



\ 



'^ RTL'^ 
V 



f = f 



V 

rT 



RT 



a(3 



RTL'^ 



where L is the characteristic length scale, r is the char- 
acteristic time scale for diffusion in the liquid, RT is 
the characteristic energy scale, and V is molar volume. 
In this work Kij and A^/s were taken to be diagonal 
and constant, and the following parameters were used: 
A = 16, Va = 1, 1^12 = ^13 = %3 = -10, Wa^ = .2, 
K = A = 8. A large negative Q insures that phases 
have limited solubility and there is a large energy barrier 
between phases in composition-space. Choosing dimen- 
sional units of L = 1 nm, V = 10 cm^, r = 1.6 x 10~^ s, 
n = 1 mol, and RT = 8.314 kJ/mol, the diffusivity in the 
liquid is D = 10~^ cm^/s, and the equilibrium interfacial 
width is about 8nm in c and 4nm in cf). The surface 
energy between solid phases is ass = 1.3 J/m^, and the 
solid-liquid surface energy is asi = .7 J/m^ at Afm = 0. 
All simulations were performed on a computational 
grid of 512 x 512 points using a time-adaptive pseudo- 
spectral method [44 that included Langevin noise. Be- 
cause the simulations are 2D, the system is effectively a 
thin film. 



If there are assumed to be no compositional effects that con- 
tribute to asymmetry in latent heat AHm, the change in free 
energy upon solidification is related to undercooling for small 
AT = Tm-T according to: 



Afrr 



AHmAT 



(35) 




(a) Below the melting point, Afm = 1.45. 




/ ^3 



C^ C2 

(b) Above the melting point, A/^ 




.2. 



FIG. 4: (Color online) Ternary eutectic free energy 

surfaces, viewed from below, and the corresponding 

phase diagrams. The diffusing components are Ci (red), 

C2 (green) and C3 (blue), and there are three solid 

phases: a red-rich phase, a green-rich phase, and a 

blue-rich phase. A silver liquid phase appears in the 

center of the phase diagram when A/^ < . 



A. Nucleation and groAvth 

Nucleation and growth in a ternary eutectic was sim- 
ulated in figure |5j The undercooling of the system cor- 
responds to Afm = 1.45 (figure 4a). The initial condi- 



tion was homogeneous metastable liquid of composition 
c = (.35, .31, .34), and all three solid phases exist at ap- 
proximately equal volume fractions at equilibrium. 

Langevin noise was added to the composition variables, 
and circular seed nuclei were added to the phase vari- 
ables. The energy of these nuclei followed a Gaussian 
distribution, and the radii was estimated using a classi- 
cal nucleation approach as described in [44]. The energy 
distribution and frequency of these nuclei was chosen so 
that phase transformation occurs in a reasonable amount 
of simulation time, and therefore is not rigorous. Given 
the size of the system being simulated, the nucleation 
rate is quite large, corresponding to a system with a high 
density of heterogeneous nucleation sites. 

Since the system is slightly enriched in ci (red compo- 
nent), the red-rich phase has the lowest nucleation barrier 
and is observed to nucleate first, as predicted by the rule 
of Stranski and Totomanow. The growing nuclei then un- 
dergo secondary nucleation at the growth front and blue 
and green solid phases are observed to form. Eutectic 
colony morphologies are expected in dilute ternary sys- 
tems [45], but there virtually no theory for multiphase 
morphology in concentrated systems, as noted in [T|[24]. 
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(a) i=15 



(b) i=30 




(c)t 



(d)t: 



FIG. 



5: (Color online) Simulation of nucleation, growth, 
and coarsening of a three-phase solid from a 
homogeneous metastable liquid. 



Figure |5] indicates that complex pattern formation is pos- 
sible when three-phase solidification is confined to a 2D 
film and all phases occur at approximately equal volume 
fractions at equilibrium. 



B. Premelting and metastable liquid 

The completely solidified structure in figure [5] was then 
brought to a higher temperature that was still below the 
melting point, and allowed to coarsen. Because isother- 
mal conditions are assumed, the temperature increase 
happens instantaneously. Figure [6] shows the system just 
before and just after the temperature increase. The mag- 
nitude of the temperature increase was not large enough 
to produce a stable liquid region in the phase diagram, 
but liquid is observed to form anyway, pooling first at 
phase triple junctions and then wetting phase bound- 
aries as temperature is increased. This behavior is qual- 
itatively similar to experimental observations [3l|4j. 

Raj [42 theorized that forming liquid at a triple junc- 
tion reduces curvature and places the liquid under neg- 
ative pressure, creating a stable melt pocket. Here we 
find that the shape and positions of the free energy sur- 
faces also plays an important role as well. The liquid 
surface lies above the convex hull of the solid curves, but 
below the solid surfaces themselves over a large compo- 
sition range. Premelting may be understood as the sys- 
tem making use of these metastable liquid states, first at 




(a) Just before the temperature increase. 




(b) Shortly after the temperature increase, i = 10. 




(c) The hquid phase fraction of 

the system at i = 10. Liquid is 

white and sohd is black. 

FIG. 6: (Color online) Premelting is observed at triple 

junctions and phase boundaries when a three-phase 

solid is heated to A/^ = -3, slightly below the melting 

point. 



triple junctions where the energy difference between the 
solid and metastable liquid surfaces is largest, and then 
at grain boundaries where the difference is smaller but 
still favors liquid over compositions far from the single 
phase solid regions. 

The composition mapqj in figure [6] reveal the effect 
of the metastable liquid surface when temperature is in- 



^ Phase-field simulations contain a lot of important information 
that must be extracted from images of microstructure. To ad- 
dress this difficulty, a composition map was developed to visually 
reveal information about compositions in a ternary microstruc- 
ture. The composition map is a triangle drawn to correspond to 
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melt when they reached a size of 10-15 nm, approximately 
twice the width of the interface. The white pockets that 



(a) Premelting is observed at 

phase triple junctions and phase 

boundaries. 



(b) Premelting is only observed 
at phase triple junctions. 



FIG. 7: (Color online) The shape of the liquid free 

energy surface affects whether premelting occurs at 

phase boundaries. The black lines denote a low energy 

path connecting the bulk compositions of two phases. 

Phase boundary premelting occurs when the black line 

traverses part of the liquid free energy surface. 



creased. In figure [6a| the composition variation at the 
diffuse solid interfaces shows up as straight, diffuse lines 
that connect the single phase regions. But when liq- 
uid forms at phase triple junctions and phase boundaries 
in figure [6bj the interfacial composition profiles bow in- 
ward toward the center of the composition map. The liq- 
uid free energy surface attracts interfacial compositions, 
and the trajectory of the interface through composition 
space changes so as to accommodate the low energy liq- 
uid states. 

The shape of the metastable regions of free energy sur- 
faces offers a possible explanation for why premelting 
is not always observed experimentally, and sometimes 
observed at triple junctions but not phase boundaries. 
When points at a diffuse interface are forced to choose be- 
tween several high energy states, the shape and position 
of free energy surfaces become important, as illustrated in 
figure [7j When the low energy path does not traverse the 
liquid surface, the system either does not premelt or must 
pay an energy penalty to adjust its trajectory to accom- 
modate liquid states. Furthermore, thickening of liquid 
film with increasing temperature, which is observed ex- 
perimentally, may be rationalized as an increasing traver- 
sal length along the liquid surface as it descends. 



C. Instability of small particles 

Coarsening theory predicts that the radius of shrink- 
ing particles will smoothly decrease to zero. However, as 



figure 6c illustrates, shrinking particles were observed to 



the phase diagram, with ci (red) at the top vertex, C2 (green) at 
the lower left, and C3 (blue) at the lower right. For every compo- 
sition in the microstructure, a corresponding point is drawn on 
the composition map. The color of each point matches the color 
of that composition in the microstructure. For clarity, composi- 
tions that are in the liquid phase are colored silver. 



are apparent in figure 6c are locations where small parti- 
cles melted. These melt pockets are temporary, and are 
eventually consumed by the surrounding solid. 

The melting of small particles is in agreement with 
analysis Wagner [46 , who found that at a given under- 
cooling, there is a critical radius below which nanocrys- 
talline materials become unstable and melt due to geo- 
metrical effects. Applying the analysis to Afm = .3 and 
our simulation parameters, we calculate the diameter of 
a critical particle surrounded by triple junctions [42 to 
be 8.4 nm, which is comparable to what was observed in 
our simulations. The discrepancy might be a result of the 
assumptions of sharp interfaces, constant surface energy, 
and an overly simple expression for A/^ in the analysis. 
A thorough understanding of the effect of diffuse inter- 
faces on the stability of multi-junctions is left for future 
work. 



D. Asymmetry from unequal difFusivity 

Another source of asymmetry between melting and so- 
lidification is that diffusion in a liquid is usually three 
to four orders of magnitude faster than in a solid. Dur- 
ing melting the phase that forms has high diffusivity, but 
during solidification the phase that forms has low dif- 
fusivity. It has been shown that the driving force for 
exchange of solute across the solid-liquid interface disap- 
pears when the diffusivity of the parent phase approaches 
zero [4TJ |47] . During solidification some of the driving 
force must be spent on trans-interface diffusion, while 
during melting all of the driving force goes into interface 
migration. When an alloy is cooled under nonequilibrium 
conditions and diffusion in the solid is limited, the com- 
position of the solid formed initially at the core of the 
solidifying structure is not the same as the composition 
at the outer edge of the structure. Due to nonequilib- 
rium solute distribution in rapidly solidified supersatu- 
rated solids, solutal melting below the melting point is 
possible. 

Figure [8a| shows a phase-field simulation of coarsening 
performed with slow diffusivity in the solid. The diffusiv- 
ity of each component was a linear function of the liquid 
phase fraction, and diffusivity in the solid regions was de- 
creased by three orders of magnitude. The structure that 
formed consists of smaller, rougher particles that are less 



equiaxed. The composition map in figure 8a reveals that 
solidification occurred at compositions outside the stable 
single phase regions. Once the system has frozen in a 
supersaturated state, evolution proceeds slowly because 
significant solid diffusion is required. 

Figures [8bHf| show the system shortly after being 
heated to A/^ = .3, which corresponds to a temperature 
below the melting point. The temperature increase ini- 
tially causes the regions of supersaturated solid to melt. 
Large pools of liquid form, but eventually the solid sur- 
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(a) t = 



(d) i = 100 




(b) i=5 



(e) i = 400 







;^*Vu' 



800 



FIG. 8: (Color online) Solutal melting and re-solidification of a rapidly solidified multiphase nanostructure held at 
A/m = .3, below the melting point. The diffusivity in solid is 1/1000 the diffusivity in the liquid. 



rounding these liquid regions grows back into the liquid. 
When the liquid has re-solidified, the composition map 
appears qualitatively similar to that in figure [6b) 



VI. CONCLUSION 

A diffuse interface model for microstructure with an 
arbitrary number of phases and components was de- 
rived from basic thermodynamic and kinetic principles. 
Interfaces were treated as thermodynamic entities and 
nonlinear diffusion equations for concentrated solutions 
were derived in accordance with the Gibbs-Duhem and 
Nernest-Einstein relations. A composition gradient en- 
ergy was included for the first time in a multiphase model 
to capture the effects of solute trapping, and an inhomo- 
geneous diffusion potential was introduced as the driving 
force for diffusion without a dilute solution approxima- 



tion. Inhomogeneous free energy for a multicomponent, 
multiphase system was obtained from a Taylor expansion 
that produced matrices of gradient energy coefficients. It 
was shown how the properties of each phase and compo- 
nent may be specified independently of the others, even 
when the phase fractions and mole fractions obey a mass 
constraint. A linear interpolation between free energy 
surfaces was used to avoid problematic pair- wise interac- 
tion of phases, and a multi-obstacle barrier was applied 
to permit arbitrary barriers between phases. 

The model is well-suited for studying phenomena 
where interfacial width is important, and captures de- 
tails of melting and solidification that have not previ- 
ously been modeled with phase-field methods. A nucle- 
ation barrier to solidification was observed, and melting 
in solids was found to start below the melting point at 
phase triple junctions and phase boundaries, where pock- 
ets of liquid and stable liquid films formed. Premelting 
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was the result of low-energy pathways through composi- 
tion space provided by metastable portions of free energy 
surfaces. Small particles were observed to be unstable to 
heating as predicted by theory, and the large difference 
in diffusion constants between solid and liquid was found 
to lead to solutal melting, common behavior in rapidly 
solidified alloys. 
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